Dependence on temperature and GC content of bubble length distributions in DNA 
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We report numerical results on the temperature dependence of the distribution of bubble lengths 
in DNA segments of various GC concentrations. Base-pair openings are described by the Peyrard- 
Bishop-Dauxois model and the corresponding thermal equilibrium distributions of bubbles are ob- 
tained through Monte Carlo calculations for bubble sizes up to the order of a hundred base-pairs. 
The dependence of the parameters of bubble length distribution on temperature and the GC content 
is investigated. We provide simple expressions which approximately describe these relations. The 
variation of the average bubble length is also presented. Finally, we find a temperature dependence 
of the exponent c in the probability distribution of bubble lengths. If an analogous dependence 
exists in the loop entropy exponent of real DNA, it may be relevant to understand overstretching 
in force-extension experiments. 
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I. INTRODUCTION 

Local openings of the DNA double helix are required 
in several biological functions, for instance transcription 
and replication. These local separations of the two DNA 
strands are mediated by a specific machinery in the cell. 
In order to deal with such complex processes, it is neces- 
sary first to understand the interactions keeping together 
the two complementary strands within a single DNA du- 
plex, as well as the properties of fluctuating DNA open- 
ings in thermal equilibrium. This fact has been realized 
long time ago and base-pairing interactions and the ther- 
mal stability of the double helix have been conventionally 
probed by increasing temperature up to the DNA denat- 
uration transition [J, Q . 

Base-pair openings (bubbles) occur in DNA due to 
thermal fluctuations even at temperatures well below the 
melting transition. It has been speculated that they may 
play a role in the recognition of specific DNA sites by 
DNA-binding proteins [3, 0, S 0]- Recent experiments 
using a novel hairpin quenching technique have been able 
to show this bubble formation in the pre-melting regime 
and characterize it quantitatively 0, S ■ By increas- 
ing temperature these bubbles grow and more bubbles 
are nucleated, thus leading to the complete separation of 
the two strands at the denaturation transition. There- 
fore statistical properties of DNA bubbles in a wide tem- 
perature regime, extending from biological temperatures 
up to the melting transition, are of particular interest. 
The purpose of this work is the investigation of the dis- 
tribution of bubble lengths and its temperature depen- 
dence for DNA sequences containing different percent- 
ages of guanine-cytosine (GC) base-pairs. The study 
is performed in the framework of the Peyrard-Bishop- 



Dauxois (PBD) model [13], from where conjectures are 
drawn about the actual distributions of bubbles in DNA. 

If our results on the temperature dependence of this 
distribution, obtained for finite bubble lengths, remain 
qualitatively valid in the asymptotic regime (for very 
large values of bubble lengths), then there would be a 
connection with the ongoing discussion about the inter- 
pretation of the DNA overstretching observed in force- 
extension experiments. In this case the possibilities of ei- 
ther a force- induced melting or the existence of a double- 
stranded elongated DNA phase (the so called S-DNA) 
have been proposed to explain the abrupt elongation of 
DNA at forces ~ 65pN ll|. This is briefly discussed in 
the concluding section. 

In a recent study we have presented the distribution of 
bubble lengths in the PBD model of DNA at 310K and 
we found that it can be described by a power-law mod- 
ified exponential [l^ . Anharmonic interactions between 
complementary bases forming base-pairs are responsible 
for the observed non-exponential distribution. The same 
form of distribution has been derived in the framework 
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of the Poland-Schcraga model 
sents a completely different theoretical approach of DNA 
denaturation than the PBD model that we use in our 
calculations. This distribution is also found for a primi- 
tive version of the PBD model, viz. the Peyrard-Bishop 
model [l^ with linear stacking interactions, but in this 
case the characteristic values of the parameters of the 
distribution are different [l7| . 

Here we examine how the bubble length distribution 
varies with temperature and present its complete depen- 
dence on both temperature and the GC fraction of the 
DNA segment. The PBD model is used for the de- 
scription of base-pair openings, where a set of continuous 
variables ?/„ represent the base-pair displacements from 
equilibrium distance and the index n labels the base-pairs 
along the DNA chain. The potential energy of the system 
consists of two parts: the on-site interaction V{yn) within 
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each base-pair and the stacking interaction U{yn, Un+i) 
between adjacent base-pairs. A Morse potential is used 
for the on-site energy, 



V{yn)^Dn{e 



a„yr. 



1) 



(1) 



where the parameters Z?„ and a„ distinguish between GC 
and AT base-pairs {Dqc = 0.075eV, aac = 6.9A-1 for 
a GC bas6-pair and Dat — 0.05eV, (Iat — 

4.2A-1 for 

an AT pair), while a nonlinear potential describes the 
stacking interaction, 

U{yn,yn+i) - y(l + pe-''('y"+«"+^))(2;„ - y„+,f , (2) 

with K = 0.025eV/A2, p = 2, and 6 = O.SSA-i. We 
use parameter values from previous works [1, H, [H, [3l ■ 
These parameters have been originally obtained from em- 
pirical fits to experimental data [l8| and are also able 
to successfully describe other experimental situations 
0, Entropic effects of DNA are described from the 
PBD Hamiltonian when studying its thermodynamics, 
leading to an entropy driven melting transition [loL [20| . 

The efhciency of the rather simple PBD model to de- 
scribe base-pair openings in DNA has led to its extensive 
use in the literaturejiji, JE [13, H 0, 0, HE [H, [H, 
MSHHHIialilllii. in particular, our choice of 
the PBD model for the study of bubble length distribu- 
tions is motivated by the success of the model {1^ in re- 
producing experimental measurements of bubble forma- 
tion [7, 8], as well as the accurate description at a quanti- 
tative level of melting curves in short DNA segments [3] ■ 
The coarse-grained description of the model allows the 
possibility to perform calculations with long sequences 
of up to tens of thousands of base pairs [s^l- Moreover, 
the popular Nearest Neighbor model [sl], which is very 
successful in describing the melting of short oligomers, 
does not reproduce well experimental data on intermedi- 
ate states for longer sequences d, H, 113, ■ Compared to 
the more conventional modified Ising type models used in 
the study of DNA melting [3,[i3], the PBD model is qual- 
itatively similar: the on-site potential of the PBD model 
is an extension of the magnetic field in the Ising type 
models, while the nonlinear stacking potential in Eq. Q 
corresponds to the interaction between neighboring spins 
and loop entropies in Ising models. We expect our results 
for the PBD model to hold qualitatively also for the Ising 
type models, since both models yield the same shape for 
the distribution of bubble sizes [I4I . 



II. RESULTS 

Considering a random DNA sequence of a given GC 
percentage, xqc, at equilibrium at temperature T, we 
calculate the distribution per base-pair of bubble lengths 
I, P{1), by counting during Monte Carlo simulations 



of the DPB model the average occurrences of openings 
(base-pair displacements) larger than a fixed threshold 
ythres at I successive base-pairs. In each simulation a ran- 
dom sequence of 1000 base pairs is used in which the AT 
or GC identity of each base pair is generated randomly 
under the constraint of a specified GC percentage. The 
properties of such random long sequences are not dif- 
ferent from natural sequences, as we showed in previous 
work 1J| comparing the results for a segment from Es- 
cherichia Colics (;aZ promoter with those from a randomly 
generated sequence with the same GC percentage. 

The Monte Carlo simulation is performed using the 
Metropolis algorithm [3^, which is used first for ther- 
malization and then for measurement runs. Results are 
averaged over several realizations using the same random 
sequence and different initial conditions for the inter-base 
distance within each base-pair and for the random num- 
ber generator. Other details of the simulation (number 
of realizations and Monte Carlo steps) are as in Ref. [13] . 
Since the used sequences are so long, the results are in- 
dependent of the precise realization of the random se- 
quence. We have run further simulations with different 
random sequences to assure this point. There is an is- 
sue about the dissociation observed in the PBD model in 
Monte Carlo simulations. In particular, due to the upper 
bound of the on-site potential for large positive displace- 
ments, for long enough simulations a complete dissoci- 
ation would eventually occur at any temperature, even 
below the melting transition 13, 331 . However, because 
of the large DNA segments considered here, the probabil- 
ity of a complete dissociation during the simulation time 
is so low that no such events were observed at any of the 
studied temperatures below the melting temperature. 

The threshold value considered for the openings is 
ythres = 1.5 A. Results for different values of this thresh- 
old are qualitatively similar. We use periodic bound- 
ary conditions in DNA segments of length 1000 base- 
pairs (sufficiently larger than the studied bubble sizes) 
and therefore our results refer to internal bubbles in long 
DNA chains. The results are independent of the length 
of the studied sequence, provided that it is longer than 
the longest observed bubbles [l^l . This has been verified 
by simulations on DNA segments of different length, re- 
sulting in the same distributions P{1) for the bubble sizes 
studied here, i.e. with I up to about a hundred base-pairs. 
For short molecules, the actual sequence does play a role, 
since it can trigger local openings in different regions at 
the end or in the middle of the sequence [3, [§] ■ In this 
case, the border effects associated can not be character- 
ized only by the GC percentage of the short sequence, 
but this problem is out of the scope of the present work. 

In Figure[l]we show bubble length distributions at var- 
ious temperatures for a GC percentage of xgc — 50% and 
87.5%. Similar plots have been obtained for nine differ- 
ent values of xgc : 0%,12.5%,25%,37.5%,...,100%. Note 
here that although the PBD model takes into account 
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Bubble length (1) 



FIG. 1: (Color online) Distribution per base-pair of bub- 
ble lengths I (in number of base-pairs), P{1), for different 
temperatures (points, as indicated in the plots) at random 
DNA sequences with a GC content of (a) 50 % and (b) 87.5 
%. Continuous lines are fits with the distribution of Eq.Q. 

ythres = I.5A. 



cooperativity in bubble formation through the nonlinear 
stacking interaction, it does not contain any parameter 
describing a bubble nucleation size. Thus the distribu- 
tions in figure [1] show results even for bubbles of length 
I — 1, implying that the nucleation size is 1. This should 
not be confused with the minimum length of a short DNA 
sequence necessary to sustain bubble states, which has 
been shown experimentally 0, [8| and theoretically within 
the DPB model [l^ to be greater than one. 

Except for the cases of i = 1 at relatively higher tem- 
peratures (closer to the melting transition), the numer- 
ical results obtained for the distributions can be rather 
well described by the power-law modified exponential [l^ 



^(0 = 1^-^, for/>l. 



(3) 



In the following we characterize the dependence of the 
parameters of distribution on T and xqc, and pro- 
vide approximate expressions for the relations ^(T, xgc)i 
c{T,xgc), and W(T,XGc)- The distribution parameters 
are obtained through fitting of plots like those of Figure 
[T]with Eq. ([3|), using a weight proportional to 1/P(/)^. 

We note here that Eq. ^ can be derived from poly- 
mer physics 1^ 14, 1^ as an asymptotic expression for 



very large bubble sizes, Z ^ 1, in DNA. However, we 
find that the same expression describes well the bubble 
length distribution in the PBD model of DNA, even for 
small bubble sizes I . We emphasize that we are interested 
here in bubble lengths up to about a hundred, or a few 
hundred at most, base-pairs (which are relevant for any 
practical purpose) and not for the asymptotic behavior 
of the distribution. In this context we use Eq. ([3]) as 
an empirical formula valid for finite bubble lengths I (in 
the non- asymptotic regime), and describe the variation 
of its parameters on T and xqc If one is interested in 
the asymptotic behavior of the distribution, this can be 
obtained by proper scaling analysis and, as it has been 
shown in Ref . (29| , it may be described by different val- 
ues of parameter c. Therefore we are not concerned with 
the order of the melting transition and the exponent c 
presented here is not indicative of the kind of the transi- 
tion 14, 15, 22, 2^, 3^, as it also depends on the some- 
how arbitrary value of the threshold chosen to consider 
a base-pair open. 

Figure [2] presents the dependence of the decay length 
^ of Eq. ([3]). In [2^ the variation of ^ with tempera- 
ture is shown (points) for different values of xqc- The 
T-dependence is accurately described by the divergent 
function 



= eo + 



(4) 



where Tc is the denaturation temperature. Such a re- 
lation is also valid for the Poland-Scheraga model [3]. 
Lines in Figure [2^ show fittings of the ^ (T) data with 
Eq. using a weight proportional to 1/^^. The criti- 
cal temperature as obtained from the fitting at differ- 
ent values of xgc is presented with circles in Figure [2Jd. 
A linear dependence of Tc on the GC content is found, 
in accordance with known experimental results (38| and 
calculations from simplified models [s^. A least square 
fitting of the Tc{xgc) data results in the continuous line 
shown in Figurej^b- Regarding the homopolymer cases of 
poly(dA)-poly(dT) {xgc = 0%) and poly(dG)-poly(dC) 
{xgc = 100%), the critical temperatures for the tran- 
sition can be independently calculated through the nu- 
merically exact transfer integral technique [io, |4l|, |43 |. 
and the corresponding results are Tc{xgc = 0) = 325. 2K 
and Tc{xgc — 100) = 366. OK. These values are shown 
with squares joined by a dashed line in Figure [2}d, where 
the line lies inside the error interval of the Monte Carlo 
results. The actual denaturation temperatures obtained 
through the Monte Carlo simulations are in agreement 
with those presented in Fig. \Bp. Therefore, the critical 
temperatures Tc shown in Fig. [2Jd provide a guide of how 
far from the melting transition are the data at various 
temperatures presented in this work. 

The other parameters ^1 and ^0 resulting from the fit- 
ting of the ^ (r) data with Eq. ^ are shown with circles 
in Figures [2t and[2Ji, respectively. Their dependence on 
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FIG. 2: (Color online) (a) Dependence of the decay length ^ 
of the distribution ([3]) on the temperature, for different val- 
ues of the GC content of the DNA sequence (points). Lines 
show fits with the function of Eq. Q. (b) Dependence of 
the critical temperature Tc, as obtained from the fitting of 
the C(r) data with Eq. (gj, on the GC content of the DNA 
sequence (circles). Solid line represents a least square fit ac- 
cording to a linear dependence. Squares show exact results 
of the critical temperatures for the homogeneous cases of 0% 
GC and 100% GC, obtained from transfer integral calcula- 
tions, while the dashed line connects these two points, (c) 
and (d) Dependence of the parameters and ^o, respectively, 
of the fit of the ^(T) data with Eq. Q, on the GC content 
of the sequence (circles). Solid lines represent linear fits of 
the corresponding data. Equations of straight lines resulting 
from the corresponding fittings are shown in (b), (c), and (d), 
where the values in parentheses represent errors of the fitting 
parameters. 
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FIG. 3: (Color online) (a) Dependence of the exponent c of the 
distribution Q on the temperature, for different values of the 
GC content of the DNA sequence (points). Lines show linear 
fits with equation (b) and (c) Dependence of the parame- 
ters ci and Co, respectively, of the fit of the c(T) data with Eq. 
((6]), on the GC content of the sequence (filled points). Error 
bars are standard errors resulting from the fitting procedure. 
Solid lines represent linear fits of the corresponding data and 
the resulting equations are also shown. 



xqc can be approximately considered as linear (contin- 
uous lines in Figs. [2]; andUJi). Therefore, the relation 



ii^^XGc) = ai + a-iXGC 



as + ajXGC 
as + aeXGC - T ' 



(5) 



where ai,a2,---,aQ are constants, can approximately 
provide the dependence of the decay length ^ on T and 

XGC- 

The variation of the exponent c of the distribution 
([3]) is presented in Figure [21 Points in [3^ show the 
temperature dependence of c for different GC percent- 
ages (for clarity of the plot the corresponding results for 
XGc = 12.5%, 37.5%, 62.5%, 87.5% have been omitted). 
This dependence may approximately be described by a 
linear function 



c(T) = Co + ciT. 



(6) 



Lines in Figure [3^ show fittings of the numerical results 
with the above formula. The parameters ci and cq ob- 
tained from the fitting at different values of GC content 
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FIG. 4: (Color online) (a) Dependence of the pre-exponential 
coefficient W of the distribution (|3]) on the temperature, 
for different values of the GC content of the DNA sequence 
(points). Lines show quadratic fits with equation Q. (b), 
(c), and (d) Dependence of the parameters W2, W\, and Wo, 
respectively, of the fit of the W{T) data with Eq. ([8]), on the 
GC content of the sequence (circles). Solid lines represent lin- 
ear fits of the corresponding data and the resulting equations 
are also shown. 



are shown with points in Figures[3jD and[5}:, respectively, 
while the corresponding error bars are derived from the 
fitting procedure. The latter plots indicate a linear de- 
pendence of ci and co on xgc^ implying the approximate 
relation 

c(T, XGc) = 61 + b2XGc + {h + biXGc)T, (7) 

where 61, 62, 63, and 64 are constants independent of T 
and xgg- 

In Figure 2] is shown the dependence of the coefficient 
W of the distribution ([3]). Points in [4^ display the vari- 
ation of W with T for various GC contents. A quadratic 
function 



WiT) = Wo + WiT + W2T^ 



(8) 



can describe rather well the temperature dependence of 
W, at least in the studied temperature regime. The cor- 
responding fittings with Eq. ([5]) are shown by lines in 
Figure [3^. The resulting fitting parameters at different 



GC percentages are plotted by points in Figures [4Jd, [41;, 
and[4ji, respectively. These plots can approximately be 
described by hnear dependences of W2, Wi, and Wq on 
xgc (solid lines). As a result the expression 

W{T,XGc) = di+d2XGc+id3+d4XGc)T+{d5+deXGc)T^ 

(9) 

can approximate the dependence of on T and xgc, 
where di , o?2 , • • ■ , ^6 are constants. 

The detailed results of this investigation do not con- 
firm a bilinear dependence of the parameters of the dis- 
tribution on the GG content at a fixed temperature [13] ■ 
Instead, a linear dependence of the exponent c and the 
coefficient W on xgc arises at constant T. We note that 
the results presented in Figure 2 of Ref. can be de- 
scribed weU by Eqs. (0, (0), and for T = 310K and 
using the values of constants as they provided in Figures 
[2jD,c,d (for ai),[3)3,c (for bi), and|3t>,c,d (for di). 

We finally present the variation of the average bubble 
length, Lb, on the two parameters of interest, viz. tem- 
perature and GG percentage. is obtained through 
the total number of base-pairs in bubble states divided 
by the total number of bubbles : 



Ei>im' 



(10) 



Substituting P{1) by the power-law modified exponential 
function of Eq. the sums in Eq. fTU]) yield: 



Lic-i(e-i/0 
Lic(e-i/«) ' 



(11) 



where Lis(z) is the first branch of the polylogarithm func- 
tion [isl, defined as Lis{z) = J2kLi 1^". The average 
bubble length for a given temperature and GG percent- 
age, Lb{T,xgc), can be obtained substituting in Eq. pT|) 
the parameters ^ and c, calculated for the corresponding 
T and xgc through Eqs. (O and respectively, using 
the values of constants determined in Figures [2] and [31 

Figure [5^ depicts the dependence oi Lb on xgc for 
different temperatures, obtained directly from our nu- 
merical simulations (points). Eq. pT|) . shown by dotted 
lines in the plot, reproduces correctly the behavior of the 
data. However, the errors in the estimated parameters 
in Figures [3 and [31 add up to produce a small deviation 
between the numerical points and the values obtained by 
Eq. pT|) . We also show, using solid fines in the plot, 
more accurate fits of the numerical data with a simpler 
phenomenological expression, namely the exponentially 
decaying function 



Lb = To + Fi exp(-r2XGc)- 



(12) 



These empirical fits generalize the apparent exponential 
decay seen earlier at 310K [13] ■ Presenting these data 
on a different way, more appropriate for experimental in- 
vestigations, the temperature dependence of Lb for fixed 
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FIG. 5: (Color online) (a) Dependence of the average bub- 
ble length L_g, Eq. (|10|l . on the GC content, for different 
temperatures (points). Dotted lines show the analytical re- 
sult of Eq. pif) . using the values of ^ and c from Eqs. ([5]) 
and (O. Solid lines show empirical fits with exponential func- 
tions, Eq. (|12p . (b) Dependence of the average bubble length 
Lb on temperature for different GC contents (points). Lines 
represent empirical fits with cubic functions. 

GC contents is shown in figure [SJd (points), for tempera- 
tures sufficiently below the melting transition. Here for 
clarity of the plot we do not show again the prediction 
of Eq. pT|) . but instead we show with lines fits with a 
simple phenomenological function, namely a cubic poly- 
nomial, which describes approximately the data in the 
investigated temperature regime. The coefficients of the 
cubic polynomials show an exponential dependence, of 
the form p^ . on the GC percentage. 



III. CONCLUSIONS 

In conclusion, we have presented the dependence of 
bubble length distributions in the PBD model of DNA 
on temperature and the GC content. The investigated 
temperature regime was extended from biologically rel- 
evant values up to values below the melting transition. 
Approximate expressions have been obtained for the pa- 
rameters of the power-law modified exponential distribu- 
tion in the non-asymptotic regime (for bubble lengths up 
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FIG. 6: Extrapolated dependence of the exponent c at the 
critical temperature (cc) on the GC content. 

to about a hundred base-pairs) . The exponent c behaves 
linearly both in temperature and GC content, Eq. ([7|. 
while the coefficient W shows a quadratic dependence on 
temperature and a linear dependence on the GC frac- 
tion, Eq. ([9]). The decay length ^ is described by the 
relatively simple equation ([5]). The constants a^, bi, and 
di appearing in these expressions depend on the ampli- 
tude Uthres of the Considered base-pair openings, and for 
Uthres = 1.5 A are given by the values shown in Figures 
m [3l and [H Using the expressions of the exponent c 
and the decay length ^, the average bubble length is an- 
alytically given by Eq. (jlip . Our results may be useful 
in biotechnological applications which involve thermally 
induced DNA denaturation, as well as in recent hairpin 
quenching experiments 0, H, @] for studying bubble for- 
mation. 

An important result of this work is the prediction of a 
dependence both in temperature and genetic sequence of 
the exponent c of the bubble distribution ([3]). This pre- 
diction can be ex per imentally verified using the method 
proposed in Ref . [4J| : based on the measurement of cer- 
tain correlation functions of base-pair openings using flu- 
orescence correlation spectroscopy (45| . the exponent c 
can be experimentally calculated. 

The sequence dependence of c is specially relevant at 
the melting temperature. At this temperature very large 
bubbles appear and the asymptotic behavior of the bub- 
ble length distribution is necessary. Assuming that the 
behavior revealed in our findings qualitatively holds in 
the asymptotic regime (which needs to be investigated 
since there is not a priori any reason for this to be true), 
then the particular dependence of c can be extrapolated. 
Combining an expression for c as that given in Eq. ([7]) and 
the dependence of critical temperature on the GC per- 
centage presented in Fig.[2}3, we show in Fig. [6] estimated 
values of c at the melting temperature, Cc, as a function 
of the GC percentage. We see that Cc is not constant, 
but decreases as the GC fraction increases. Although, 
in view of the above assumptions, this result should be 
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taken with caution, such a behavior is consistent with the 
observation, coming from lattice models simulations that 
consider DNA strands as self- avoiding random walks, of 
a decrease of Cc with increasing stiffness of the strands 

Under the same assumption that our results for the 
temperature dependence of c qualitatively hold in the 
asymptotic regime, then this may be useful to un- 
derstand experimental results in force-extension exper- 
iments. When performing single molecule experiments, 
force-extension curves of double stranded DNA display 



hysteresis [47|, |48| . The hysteresis observed is more pro- 
nounced at higher temperatures. The existence of an 
elongated double-stranded form of DNA, called S-DNA, 
has been proposed to explain the abrupt elongation of 
DNA at a force of ^ 65vN [491 . Numerical studies of the 
force-extension problem [SO'I suggest that the existence 
of S-DNA is necessary to explain the temperature depen- 
dence of the observed hysteresis. However, an alternative 
explanation without the appearance of an S-DNA phase 
possible if the exponent c is temperature dependent 
The order of variation of c we find here is in the 
right direction (c increases with temperature, implying 
larger hysteresis loops at high temperatures according to 
the experiments) and of the necessary magnitude [50l | to 
explain the hysteresis effects observed. Hence, if the qual- 
itative behavior of our results is valid in the asymptotic 
regime, this would have an effect on the ongoing debate 
about S-DNA, favoring force-induced melting as an ex- 
planation for the DNA elongation over the formation of 
an S-DNA phase. Interestingly, recent experimental re- 
sults support the force-induced melting explanation 5l( . 
in agreement with the expectation from the temperature 
dependence of c. 
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